using StructuralIdentifiability
ode =  @ODEmodel(
                   x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
                   x2'(t) = alpha * x1(t) - beta * x2(t),
                   x3'(t) = gama * x2(t) - delta * x3(t),
                   x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
                   y(t) = x1(t)
               )
find_identifiable_functions(ode, with_states=true)
#=
 ...
 (alpha*gama)//x3(t)
 (x2(t)*gama - x3(t)*delta)//x3(t)
=#

#=
 ...
 x3(t)//(alpha*gama)
 (x2(t)*gama - x3(t)*delta)//x3(t)
=#